
#Replication of all Empirical Analyses in the manuscript and the Appendix
#Full data drawn from other sources is includes to make extensions of the replication easier.
#V-Dem, PTS, Polity, and Shen & Truex's Replication code are included among the data files

rm(list = ls())

library(countrycode)
library(ggplot2)
library(ggridges)
library(ggjoy)
library(ggpubr)
library(lme4)
library(ggeffects)
library(GGally)
library(dotwhisker)
library(survival)
library(Zelig)
library(foreign)
library(purrr)
library(stringr)
library(Rcpp)
library(Hmisc)
library(dplyr)
library(estimatr)
library(effects)
library(redres)
library(texreg)
library(stargazer)

#Data Recoding 1
#Code is inefficient, but was easier to track in relation to the codebook and check for any errors.

wvs <- readRDS("WVSo.rds")
wvs <- data.frame(wvs)

#Coding which data should be included (removing nonresponses associated with questions not being asked, etc.)

wvs$govOB <- NA
wvs$govOB[wvs$E069_11==-1] <- 1
wvs$govOB[wvs$E069_11==-2] <- 1
wvs$govOB[wvs$E069_11==-4] <- NA
wvs$govOB[wvs$E069_11==-5] <- NA
wvs$govOB[wvs$E069_11==1] <- 1
wvs$govOB[wvs$E069_11==2] <- 1
wvs$govOB[wvs$E069_11==3] <- 1
wvs$govOB[wvs$E069_11==4] <- 1

wvs$hrOB <- NA
wvs$hrOB[wvs$E124==-1] <- 1
wvs$hrOB[wvs$E124==-2] <- 1
wvs$hrOB[wvs$E124==-3] <- NA
wvs$hrOB[wvs$E124==-4] <- NA
wvs$hrOB[wvs$E124==-5] <- NA
wvs$hrOB[wvs$E124==1] <- 1
wvs$hrOB[wvs$E124==2] <- 1
wvs$hrOB[wvs$E124==3] <- 1
wvs$hrOB[wvs$E124==4] <- 1

wvs$demOB <- NA
wvs$demOB[wvs$E236==-1] <- 1
wvs$demOB[wvs$E236==-2] <- 1
wvs$demOB[wvs$E236==-3] <- NA
wvs$demOB[wvs$E236==-4] <- NA
wvs$demOB[wvs$E236==-5] <- NA
wvs$demOB[wvs$E236==0] <- 1
wvs$demOB[wvs$E236==1] <- 1
wvs$demOB[wvs$E236==2] <- 1
wvs$demOB[wvs$E236==3] <- 1
wvs$demOB[wvs$E236==4] <- 1
wvs$demOB[wvs$E236==5] <- 1
wvs$demOB[wvs$E236==6] <- 1
wvs$demOB[wvs$E236==7] <- 1
wvs$demOB[wvs$E236==8] <- 1
wvs$demOB[wvs$E236==9] <- 1
wvs$demOB[wvs$E236==10] <- 1

wvs$lifeOB <- NA
wvs$lifeOB[wvs$A170==-1] <- 1
wvs$lifeOB[wvs$A170==-2] <- 1
wvs$lifeOB[wvs$A170==-3] <- NA
wvs$lifeOB[wvs$A170==-4] <- NA
wvs$lifeOB[wvs$A170==-5] <- NA
wvs$lifeOB[wvs$A170==0] <- 1
wvs$lifeOB[wvs$A170==1] <- 1
wvs$lifeOB[wvs$A170==2] <- 1
wvs$lifeOB[wvs$A170==3] <- 1
wvs$lifeOB[wvs$A170==4] <- 1
wvs$lifeOB[wvs$A170==5] <- 1
wvs$lifeOB[wvs$A170==6] <- 1
wvs$lifeOB[wvs$A170==7] <- 1
wvs$lifeOB[wvs$A170==8] <- 1
wvs$lifeOB[wvs$A170==9] <- 1
wvs$lifeOB[wvs$A170==10] <- 1

wvs$trustOB <- NA
wvs$trustOB[wvs$A165==-1] <- 1
wvs$trustOB[wvs$A165==-2] <- 1
wvs$trustOB[wvs$A165==-4] <- NA
wvs$trustOB[wvs$A165==-5] <- NA
wvs$trustOB[wvs$A165==1] <- 1
wvs$trustOB[wvs$A165==2] <- 1

wvs$tvOB <- NA
wvs$tvOB[wvs$E069_10==-1] <- 1
wvs$tvOB[wvs$E069_10==-2] <- 1
wvs$tvOB[wvs$E069_10==-4] <- NA
wvs$tvOB[wvs$E069_10==-5] <- NA
wvs$tvOB[wvs$E069_10==1] <- 1
wvs$tvOB[wvs$E069_10==2] <- 1
wvs$tvOB[wvs$E069_10==3] <- 1
wvs$tvOB[wvs$E069_10==4] <- 1



#Do not know and refuse to answer marked as sensitive.

wvs$govDK <- NA
wvs$govDK[wvs$E069_11==-1] <- 1
wvs$govDK[wvs$E069_11==-2] <- 1
wvs$govDK[wvs$E069_11==-4] <- NA
wvs$govDK[wvs$E069_11==-5] <- NA
wvs$govDK[wvs$E069_11==1] <- 0
wvs$govDK[wvs$E069_11==2] <- 0
wvs$govDK[wvs$E069_11==3] <- 0
wvs$govDK[wvs$E069_11==4] <- 0

wvs$hrDK <- NA
wvs$hrDK[wvs$E124==-1] <- 1
wvs$hrDK[wvs$E124==-2] <- 1
wvs$hrDK[wvs$E124==-3] <- NA
wvs$hrDK[wvs$E124==-4] <- NA
wvs$hrDK[wvs$E124==-5] <- NA
wvs$hrDK[wvs$E124==1] <- 0
wvs$hrDK[wvs$E124==2] <- 0
wvs$hrDK[wvs$E124==3] <- 0
wvs$hrDK[wvs$E124==4] <- 0

wvs$demDK <- NA
wvs$demDK[wvs$E236==-1] <- 1
wvs$demDK[wvs$E236==-2] <- 1
wvs$demDK[wvs$E236==-3] <- NA
wvs$demDK[wvs$E236==-4] <- NA
wvs$demDK[wvs$E236==-5] <- NA
wvs$demDK[wvs$E236==0] <- 0
wvs$demDK[wvs$E236==1] <- 0
wvs$demDK[wvs$E236==2] <- 0
wvs$demDK[wvs$E236==3] <- 0
wvs$demDK[wvs$E236==4] <- 0
wvs$demDK[wvs$E236==5] <- 0
wvs$demDK[wvs$E236==6] <- 0
wvs$demDK[wvs$E236==7] <- 0
wvs$demDK[wvs$E236==8] <- 0
wvs$demDK[wvs$E236==9] <- 0
wvs$demDK[wvs$E236==10] <- 0

#Preserve Actual Values
#Missing as NA
wvs$govMissNA <- wvs$E069_11
wvs$govMissNA[wvs$E069_11==-1] <- NA
wvs$govMissNA[wvs$E069_11==-2] <- NA
wvs$govMissNA[wvs$E069_11==-4] <- NA
wvs$govMissNA[wvs$E069_11==-5] <- NA
wvs$govMissNA <- 4 - wvs$govMissNA

wvs$hrMissNA <- wvs$E124
wvs$hrMissNA[wvs$E124==1] <- 0
wvs$hrMissNA[wvs$E124==2] <- 1
wvs$hrMissNA[wvs$E124==-1] <- NA
wvs$hrMissNA[wvs$E124==-2] <- NA
wvs$hrMissNA[wvs$E124==-3] <- NA
wvs$hrMissNA[wvs$E124==-4] <- NA
wvs$hrMissNA[wvs$E124==-5] <- NA
wvs$hrMissNA <- 4 - wvs$hrMissNA

wvs$demMissNA <- wvs$E236
wvs$demMissNA[wvs$E236==-1] <- NA
wvs$demMissNA[wvs$E236==-2] <- NA
wvs$demMissNA[wvs$E236==-3] <- NA
wvs$demMissNA[wvs$E236==-4] <- NA
wvs$demMissNA[wvs$E236==-5] <- NA


#Recoding

wvs$E069_11[wvs$E069_11==1] <- 0
wvs$E069_11[wvs$E069_11==2] <- 1
wvs$E069_11[wvs$E069_11==-1] <- 2
wvs$E069_11[wvs$E069_11==-2] <- 2
wvs$E069_11[wvs$E069_11==-4] <- NA
wvs$E069_11[wvs$E069_11==-5] <- NA
wvs$E069_11 <- 4 - wvs$E069_11


wvs$E124[wvs$E124==1] <- 0
wvs$E124[wvs$E124==2] <- 1
wvs$E124[wvs$E124==-1] <- 2
wvs$E124[wvs$E124==-2] <- 2
wvs$E124[wvs$E124==-3] <- NA
wvs$E124[wvs$E124==-4] <- NA
wvs$E124[wvs$E124==-5] <- NA
wvs$E124 <- 4 - wvs$E124

wvs$E236[wvs$E236==-1] <- 5
wvs$E236[wvs$E236==-2] <- 5
wvs$E236[wvs$E236==-3] <- NA
wvs$E236[wvs$E236==-4] <- NA
wvs$E236[wvs$E236==-5] <- NA

wvs$C006[wvs$C006==-1] <- NA
wvs$C006[wvs$C006==-2]<-NA
wvs$C006[wvs$C006==-3]<-NA
wvs$C006[wvs$C006==-4]<-NA
wvs$C006[wvs$C006==-5]<-NA

#Missing Data Not Sensitive

wvs$lifeDK <- NA
wvs$lifeDK[wvs$A170==-1] <- 1
wvs$lifeDK[wvs$A170==-2] <- 1
wvs$lifeDK[wvs$A170==-3] <- NA
wvs$lifeDK[wvs$A170==-4] <- NA
wvs$lifeDK[wvs$A170==-5] <- NA
wvs$lifeDK[wvs$A170==0] <- 0
wvs$lifeDK[wvs$A170==1] <- 0
wvs$lifeDK[wvs$A170==2] <- 0
wvs$lifeDK[wvs$A170==3] <- 0
wvs$lifeDK[wvs$A170==4] <- 0
wvs$lifeDK[wvs$A170==5] <- 0
wvs$lifeDK[wvs$A170==6] <- 0
wvs$lifeDK[wvs$A170==7] <- 0
wvs$lifeDK[wvs$A170==8] <- 0
wvs$lifeDK[wvs$A170==9] <- 0
wvs$lifeDK[wvs$A170==10] <- 0

wvs$trustDK <- NA
wvs$trustDK[wvs$A165==-1] <- 1
wvs$trustDK[wvs$A165==-2] <- 1
wvs$trustDK[wvs$A165==-4] <- NA
wvs$trustDK[wvs$A165==-5] <- NA
wvs$trustDK[wvs$A165==1] <- 0
wvs$trustDK[wvs$A165==2] <- 0

wvs$tvDK <- NA
wvs$tvDK[wvs$E069_10==-1] <- 1
wvs$tvDK[wvs$E069_10==-2] <- 1
wvs$tvDK[wvs$E069_10==-4] <- NA
wvs$tvDK[wvs$E069_10==-5] <- NA
wvs$tvDK[wvs$E069_10==1] <- 0
wvs$tvDK[wvs$E069_10==2] <- 0
wvs$tvDK[wvs$E069_10==3] <- 0
wvs$tvDK[wvs$E069_10==4] <- 0




#Putting together the first dataset


wvs$n <- 1
wvs$COW <- countrycode(wvs$S003A, origin='iso3n', destination='cown', nomatch=NA)
wvs$ccode <- wvs$COW
wvs$year <- wvs$S020

wvs <- wvs %>%
  group_by(COW,S020)  %>% 
  mutate(weightedGov = weighted.mean(E069_11, S017,na.rm=T),
         weightedHR = weighted.mean(E124,S017,na.rm=T),
         weightedDem = weighted.mean(E236, S017,na.rm=T),
         weightedGovNA = weighted.mean(govMissNA, S017,na.rm=T),
         weightedHRNA = weighted.mean(hrMissNA, S017,na.rm=T),
         weightedDemNA = weighted.mean(demMissNA, S017,na.rm=T),
         weightedEconomic = weighted.mean(C006,S017,na.rm=T))

SCI <- aggregate(list("govDK"=wvs$govDK,"hrDK"=wvs$hrDK,"demDK"=wvs$demDK,
                      "tvDK"=wvs$tvDK, "lifeDK"=wvs$lifeDK, "trustDK"=wvs$trustDK,
                      "govOB"=wvs$govOB,"hrOB"=wvs$hrOB,"demOB"=wvs$demOB,
                      "tvOB"=wvs$tvOB,"lifeOB"=wvs$lifeOB,"trustOB"=wvs$trustOB),
                 by=list("ccode"=wvs$COW,"year"=wvs$S020),sum,
                 na.rm=T, na.omit=NULL)

GovSupport <- aggregate(list("Government"=wvs$weightedGov,"HR"=wvs$weightedHR,"Democratic"=wvs$weightedDem, 
                             "govDK"=wvs$govDK, "hrDK"=wvs$hrDK,"demDK"=wvs$demDK,
                             "economicCondition"=wvs$weightedEconomic,
                             "govNonNA"=wvs$weightedGovNA, "hrNonNA"=wvs$weightedHRNA,"demNonNA"=wvs$weightedDemNA),
                        by=list("ccode"=wvs$COW,"year"=wvs$S020),mean,
                        na.rm=T, na.omit=F)


GovSupport$Government[is.nan(GovSupport$Government)] <- NA
GovSupport$HR[is.nan(GovSupport$HR)] <- NA
GovSupport$Democratic[is.nan(GovSupport$Democratic)] <- NA
GovSupport$economicCondition[is.nan(GovSupport$economicCondition)] <- NA
GovSupport$govNonNA[is.nan(GovSupport$govNonNA)] <- NA
GovSupport$hrNonNA[is.nan(GovSupport$hrNonNA)] <- NA
GovSupport$demNonsNA[is.nan(GovSupport$demNonNA)] <- NA

SCI$RAzero <- SCI$govOB*SCI$hrOB*SCI$demOB*SCI$tvOB*SCI$lifeOB*SCI$trustOB
SCI$RAzero[SCI$RAzero !=0] <- 1
SCI$NSzero <- SCI$tvOB+SCI$lifeOB+SCI$trustOB
SCI$RA3zero <- SCI$govOB+SCI$hrOB+SCI$demOB
SCI$RA3zero[SCI$RA3zero !=0] <- 1
SCI$SCI <- (SCI$govDK + SCI$hrDK + SCI$demDK)/(SCI$govOB+SCI$hrOB+SCI$demOB) - 
  (SCI$tvDK+SCI$lifeDK+SCI$trustDK)/(SCI$tvOB+SCI$lifeOB+SCI$trustOB)

SCI$SCI[is.nan(SCI$SCI)] <- NA
SCI$SCIallPos <- SCI$SCI
SCI$SCI[SCI$RAzero==0] <- NA




SCI <- data.frame("year"=SCI$year,"ccode"=SCI$ccode,"SCI"=SCI$SCI,"SCIallPos"=SCI$SCIallPos)
GovSupport <- merge(GovSupport,SCI, by = c("ccode","year"), all.x=T)

polity <- read.csv("polity.csv")

polity$polity[polity$polity==-88] <- NA
polity$polity[polity$polity==-77] <- NA
polity$polity[polity$polity==-66] <- NA

load("PTS.RData")

vdem <- readRDS("VDEM.rds")
altdem <- data.frame("ccode"=vdem$COWcode,"year"=vdem$year, "FoE"=vdem$v2x_freexp_altinf)

PTS_2020$ccode <- PTS_2020$COW_Code_N
PTS_2020$year <- PTS_2020$Year

shen <- read.csv("SHENTRUEX.csv")
shen$n <- NULL
shen$country.name <- NULL
shen$polity <- NULL
shen$ccode <- shen$cown
shen$cown <- NULL


fullData <- merge(GovSupport, polity, by=c("ccode","year"), all.x = T)
fullData <- merge(fullData, PTS_2020, by=c("ccode","year"), all.x = T)
fullData <- merge(fullData, altdem, by=c("ccode","year"), all.x = T)
fullData <- merge(fullData, shen, by=c("ccode","year"), all.x = T)



#Creating different datasets with subsets based on thresholds of how much freedom of expression is allowed.

reducedData5 <- fullData[fullData$FoE > 0.5,]
reducedData6 <- fullData[fullData$FoE > 0.6,]
reducedData7 <- fullData[fullData$FoE > 0.7,]
reducedData8 <- fullData[fullData$FoE > 0.8,]
reducedData9 <- fullData[fullData$FoE > 0.9,]
medFoE <- median(fullData$FoE,na.rm = T)
reducedDataMedian <- fullData[fullData$FoE > medFoE,]
reducedDataMedianBelow <- fullData[fullData$FoE < medFoE,]
autoData <- fullData[fullData$authoritarian==1,]

#Using full data with NA as middle group

oneMod1 <- lm_robust(scale(Government) ~ scale(SCI) + scale(FoE) + scale(polity2) + 
               scale(gdppc) + scale(economicCondition) + scale(edu) + oil, 
             data = fullData)
twoMod1 <- lm_robust(scale(HR) ~ scale(SCI) + scale(FoE) + scale(polity2) + 
                      scale(gdppc) + scale(economicCondition) + scale(edu) + oil, 
                    data = fullData)
threeMod1 <- lm_robust(scale(Democratic) ~ scale(SCI) + scale(FoE) + scale(polity2) + 
                      scale(gdppc) + scale(economicCondition) + scale(edu) + oil, 
                    data = fullData)

oneMod2 <- lm_robust(scale(Government) ~ scale(SCI) + scale(FoE) + scale(polity2) + 
                      scale(gdppc) + scale(economicCondition) + scale(edu) + oil, 
                    data = reducedData5)
twoMod2 <- lm_robust(scale(HR) ~ scale(SCI) + scale(FoE) + scale(polity2) + 
                      scale(gdppc) + scale(economicCondition) + scale(edu) + oil, 
                    data = reducedData5)
threeMod2 <- lm_robust(scale(Democratic) ~ scale(SCI) + scale(FoE) + scale(polity2) + 
                        scale(gdppc) + scale(economicCondition) + scale(edu) + oil, 
                      data = reducedData5)

oneMod3 <- lm_robust(scale(Government) ~ scale(SCI) + scale(FoE) + scale(polity2) + 
                      scale(gdppc) + scale(economicCondition) + scale(edu) + oil, 
                    data = reducedDataMedian)
twoMod3 <- lm_robust(scale(HR) ~ scale(SCI) + scale(FoE) + scale(polity2) + 
                      scale(gdppc) + scale(economicCondition) + scale(edu) + oil, 
                    data = reducedDataMedian)
threeMod3 <- lm_robust(scale(Democratic) ~ scale(SCI) + scale(FoE) + scale(polity2) + 
                        scale(gdppc) + scale(economicCondition) + scale(edu) + oil, 
                      data = reducedDataMedian)

oneMod4 <- lm_robust(scale(Government) ~ scale(SCI) + 
                      scale(gdppc) + scale(economicCondition) + scale(edu) + oil, 
                    data = fullData)
twoMod4 <- lm_robust(scale(HR) ~ scale(SCI) + 
                      scale(gdppc) + scale(economicCondition) + scale(edu) + oil, 
                    data = fullData)
threeMod4 <- lm_robust(scale(Democratic) ~ scale(SCI) + 
                        scale(gdppc) + scale(economicCondition) + scale(edu) + oil, 
                      data = fullData)


#Table A5 (setting nonresponse as the middle category)

texreg(list(oneMod1,twoMod1,threeMod1,
            oneMod2,twoMod2,threeMod2,
            oneMod3,twoMod3,threeMod3,
            oneMod4,twoMod4,threeMod4),
       custom.coef.names = c("Intercept","SCI","Freedom of Expression","Polity","GDPPC", "Financial Situation", 
                             "Education","Oil"),
       custom.model.names  = c("Gov","HR","Dem","Gov","HR","Dem",
                               "Gov","HR","Dem","Gov","HR","Dem"),
       stars = c(0.05,0.01,0.001),
       reorder.coef = c(2,3,4,5,6,7,8,1),
       fontsize = "footnotesize",
       include.rsquared = F,
       include.fstatistic = F,
       include.adjrs = F,
       include.df=F,
       include.pvalues=F,
       custom.header = list("FullSample"= 1:3,"FoE $>$ 0.5" = 4:6, "FoE $>$ Median" = 7:9, "w/o FoE and Polity" = 10:12),
       file = "TableA5.tex")


#Figure 7

FoEPlot1 <- ggplot(fullData, aes(FoE, govNonNA)) + 
  geom_smooth() + geom_point() + ylab("Government Support") + xlab("Freedom of Expression") +
  theme_bw() + ggtitle("A") +
  theme(plot.title = element_text(hjust = 0.5),axis.line = element_line(colour = "black"),
        panel.border = element_blank(),
        panel.background = element_blank())
FoEPlot1

MissingPlotFoE <- ggplot(fullData, aes(FoE, SCI)) + 
  geom_smooth() + geom_point() + ylab("SCI") + xlab("Freedom of Expression") +
  theme_bw() + ggtitle("C") +
  theme(plot.title = element_text(hjust = 0.5),axis.line = element_line(colour = "black"),
        panel.border = element_blank(),
        panel.background = element_blank())
MissingPlotFoE


MissingPlotGov1 <- ggplot(fullData, aes(SCI, govNonNA)) + 
  geom_smooth() + geom_point() + ylab("Government Support") + xlab("SCI") +
  theme_bw() + ggtitle("C") +
  theme(plot.title = element_text(hjust = 0.5),axis.line = element_line(colour = "black"),
        panel.border = element_blank(),
        panel.background = element_blank())
MissingPlotGov1

MissingPlotGovDk1 <- ggplot(fullData, aes(govDK, govNonNA)) + 
  geom_smooth() + geom_point() + ylab("Government Support") + xlab("Question Nonresponse Rate") +
  theme_bw() + ggtitle("D") +
  theme(plot.title = element_text(hjust = 0.5),axis.line = element_line(colour = "black"),
        panel.border = element_blank(),
        panel.background = element_blank())
MissingPlotGovDk1


wvsLowess1 <- ggarrange(FoEPlot1,MissingPlotFoE,MissingPlotGov1,MissingPlotGovDk1)  
wvsLowess1 <- annotate_figure(wvsLowess1,
                              top = text_grob("Government Support, FoE and Missing Responses (Nonresponse as NA)", 
                                              color = "black"))
ggexport(wvsLowess1, filename = "Figure7.pdf")


#Only authoritarian states
oneMod1a <- lm_robust(scale(govNonNA) ~ scale(SCIallPos) + scale(FoE) + scale(polity2) + 
                       scale(gdppc) + scale(economicCondition) + scale(edu) + oil, 
                     data = autoData)
twoMod1a <- lm_robust(scale(hrNonNA) ~ scale(SCIallPos) + scale(FoE) + scale(polity2) + 
                       scale(gdppc) + scale(economicCondition) + scale(edu) + oil, 
                     data = autoData)
threeMod1a <- lm_robust(scale(demNonNA) ~ scale(SCIallPos) + scale(FoE) + scale(polity2) + 
                         scale(gdppc) + scale(economicCondition) + scale(edu) + oil, 
                       data = autoData)

plot1a <- dwplot(list(oneMod1a,twoMod1a,threeMod1a),vline = geom_vline(xintercept = 0, colour = "grey60", linetype = 2))  %>%
  relabel_predictors(c("scale(SCIallPos)"="SCI","scale(FoE)"="FoE Index","scale(polity2)"="Polity",
                       "scale(gdppc)"="GDPPC","scale(economicCondition)"="Economic Condition",
                       "scale(edu)"="Education","oil"="Oil")) +
  theme_bw() + ggtitle("Full Sample") + 
  theme(plot.title = element_text(hjust = 0.5),legend.background = element_rect(colour="grey80")) +
  scale_colour_grey(name = "Regime Assessment", 
                    labels = c("Government", "Human Rights","Democracy")) +
  scale_shape_discrete(name = "Model",
                       labels = c("Government", "Human Rights","Democracy"))


oneMod4a <- lm_robust(scale(govNonNA) ~ scale(SCIallPos) + 
                       scale(gdppc) + scale(economicCondition) + scale(edu) + oil, 
                     data = autoData)
twoMod4a <- lm_robust(scale(hrNonNA) ~ scale(SCIallPos) + 
                       scale(gdppc) + scale(economicCondition) + scale(edu) + oil, 
                     data = autoData)
threeMod4a <- lm_robust(scale(demNonNA) ~ scale(SCIallPos) + 
                         scale(gdppc) + scale(economicCondition) + scale(edu) + oil, 
                       data = autoData)

plot4a <- dwplot(list(oneMod4a,twoMod4a,threeMod4a),vline = geom_vline(xintercept = 0, colour = "grey60", linetype = 2)) %>%
  relabel_predictors(c("scale(SCIallPos)"="SCI","scale(gdppc)"="GDPPC","scale(economicCondition)"="Economic Condition",
                       "scale(edu)"="Education","oil"="Oil")) +
  theme_bw() + ggtitle("Without FoE and Polity") + 
  theme(plot.title = element_text(hjust = 0.5),legend.background = element_rect(colour="grey80")) +
  scale_colour_grey(name = "Regime Assessment", 
                    labels = c("Government", "Human Rights","Democracy")) +
  scale_shape_discrete(name = "Model",
                       labels = c("Government", "Human Rights","Democracy"))

lowPolity <- fullData[fullData$polity2 < 5.5,]
oneMod1p <- lm_robust(scale(govNonNA) ~ scale(SCIallPos) + scale(FoE) + scale(polity2) + 
                        scale(gdppc) + scale(economicCondition) + scale(edu) + oil, 
                      data = lowPolity)
twoMod1p <- lm_robust(scale(hrNonNA) ~ scale(SCIallPos) + scale(FoE) + scale(polity2) + 
                        scale(gdppc) + scale(economicCondition) + scale(edu) + oil, 
                      data = lowPolity)
threeMod1p <- lm_robust(scale(demNonNA) ~ scale(SCIallPos) + scale(FoE) + scale(polity2) + 
                          scale(gdppc) + scale(economicCondition) + scale(edu) + oil, 
                        data = lowPolity)

plot1p <- dwplot(list(oneMod1a,twoMod1a,threeMod1a),vline = geom_vline(xintercept = 0, colour = "grey60", linetype = 2))  %>%
  relabel_predictors(c("scale(SCIallPos)"="SCI","scale(FoE)"="FoE Index","scale(polity2)"="Polity",
                       "scale(gdppc)"="GDPPC","scale(economicCondition)"="Economic Condition",
                       "scale(edu)"="Education","oil"="Oil")) +
  theme_bw() + ggtitle("Polity < 6") + 
  theme(plot.title = element_text(hjust = 0.5),legend.background = element_rect(colour="grey80")) +
  scale_colour_grey(name = "Regime Assessment", 
                    labels = c("Government", "Human Rights","Democracy")) +
  scale_shape_discrete(name = "Model",
                       labels = c("Government", "Human Rights","Democracy"))


#Figure A8

wvsRegA <- ggarrange(plot1a,plot4a,plot1p, common.legend = T, legend = "bottom")
ggexport(wvsRegA, filename = "FigureA8.pdf")

#Figure A7

FoEPlot1a <- ggplot(lowPolity, aes(FoE, govNonNA)) + 
  geom_smooth() + geom_point() + ylab("Government Support") + xlab("Freedom of Expression") +
  theme_bw() + ggtitle("A") +
  theme(plot.title = element_text(hjust = 0.5),axis.line = element_line(colour = "black"),
        panel.border = element_blank(),
        panel.background = element_blank())
FoEPlot1a


MissingPlotGov1a <- ggplot(lowPolity, aes(SCIallPos, govNonNA)) + 
  geom_smooth() + geom_point() + ylab("Government Support") + xlab("SCI") +
  theme_bw() + ggtitle("B") +
  theme(plot.title = element_text(hjust = 0.5),axis.line = element_line(colour = "black"),
        panel.border = element_blank(),
        panel.background = element_blank())
MissingPlotGov1a

MissingPlotGovDk1a <- ggplot(lowPolity, aes(govDK, govNonNA)) + 
  geom_smooth() + geom_point() + ylab("Government Support") + xlab("Question Nonresponse Rate") +
  theme_bw() + ggtitle("C") +
  theme(plot.title = element_text(hjust = 0.5),axis.line = element_line(colour = "black"),
        panel.border = element_blank(),
        panel.background = element_blank())
MissingPlotGovDk1a

autoBiv <- ggarrange(FoEPlot1a,MissingPlotGov1a,MissingPlotGovDk1a)
ggexport(autoBiv,filename = "FigureA7.pdf")


#Main Analysis (OLS country-wave level)
oneMod8 <- lm_robust(scale(govNonNA) ~ scale(SCI) + scale(FoE) + scale(polity2) + 
                      scale(gdppc) + scale(economicCondition) + scale(edu) + oil, 
                    data = fullData)
twoMod8 <- lm_robust(scale(hrNonNA) ~ scale(SCI) + scale(FoE) + scale(polity2) + 
                      scale(gdppc) + scale(economicCondition) + scale(edu) + oil, 
                    data = fullData)
threeMod8 <- lm_robust(scale(demNonNA) ~ scale(SCI) + scale(FoE) + scale(polity2) + 
                        scale(gdppc) + scale(economicCondition) + scale(edu) + oil, 
                      data = fullData)

plot1NA <- dwplot(list(oneMod8,twoMod8,threeMod8),vline = geom_vline(xintercept = 0, colour = "grey60", linetype = 2))  %>%
  relabel_predictors(c("scale(SCI)"="SCI","scale(FoE)"="FoE Index","scale(polity2)"="Polity",
                       "scale(gdppc)"="GDPPC","scale(economicCondition)"="Economic Condition",
                       "scale(edu)"="Education","oil"="Oil")) +
  theme_bw() + ggtitle("Full Sample") + 
  theme(plot.title = element_text(hjust = 0.5),legend.background = element_rect(colour="grey80")) +
  scale_colour_grey(name = "Regime Assessment", 
                    labels = c("Government", "Human Rights","Democracy")) +
  scale_shape_discrete(name = "Model",
                       labels = c("Government", "Human Rights","Democracy"))

oneMod5 <- lm_robust(scale(govNonNA) ~ scale(SCI) + scale(FoE) + scale(polity2) + 
                      scale(gdppc) + scale(economicCondition) + scale(edu) + oil, 
                    data = reducedData5)
twoMod5 <- lm_robust(scale(hrNonNA) ~ scale(SCI) + scale(FoE) + scale(polity2) + 
                      scale(gdppc) + scale(economicCondition) + scale(edu) + oil, 
                    data = reducedData5)
threeMod5 <- lm_robust(scale(demNonNA) ~ scale(SCI) + scale(FoE) + scale(polity2) + 
                        scale(gdppc) + scale(economicCondition) + scale(edu) + oil, 
                      data = reducedData5)

plot2NA <- dwplot(list(oneMod5,twoMod5,threeMod5),vline = geom_vline(xintercept = 0, colour = "grey60", linetype = 2))  %>%
  relabel_predictors(c("scale(SCI)"="SCI","scale(FoE)"="FoE Index","scale(polity2)"="Polity",
                       "scale(gdppc)"="GDPPC","scale(economicCondition)"="Economic Condition",
                       "scale(edu)"="Education","oil"="Oil")) +
  theme_bw() + ggtitle("FoE Index > 0.5") +
  theme(plot.title = element_text(hjust = 0.5),legend.background = element_rect(colour="grey80")) +
  scale_colour_grey(name = "Regime Assessment", 
                    labels = c("Government", "Human Rights","Democracy")) +
  scale_shape_discrete(name = "Model",
                       labels = c("Government", "Human Rights","Democracy"))

oneMod6 <- lm_robust(scale(govNonNA) ~ scale(SCI) + scale(FoE) + scale(polity2) + 
                      scale(gdppc) + scale(economicCondition) + scale(edu) + oil, 
                    data = reducedDataMedian)
twoMod6 <- lm_robust(scale(hrNonNA) ~ scale(SCI) + scale(FoE) + scale(polity2) + 
                      scale(gdppc) + scale(economicCondition) + scale(edu) + oil, 
                    data = reducedDataMedian)
threeMod6 <- lm_robust(scale(demNonNA) ~ scale(SCI) + scale(FoE) + scale(polity2) + 
                        scale(gdppc) + scale(economicCondition) + scale(edu) + oil, 
                      data = reducedDataMedian)

plot3NA <- dwplot(list(oneMod6,twoMod6,threeMod6),vline = geom_vline(xintercept = 0, colour = "grey60", linetype = 2))  %>%
  relabel_predictors(c("scale(SCI)"="SCI","scale(FoE)"="FoE Index","scale(polity2)"="Polity",
                       "scale(gdppc)"="GDPPC","scale(economicCondition)"="Economic Condition",
                       "scale(edu)"="Education","oil"="Oil")) +
  theme_bw() + ggtitle("FoE Index > Median") +
  theme(plot.title = element_text(hjust = 0.5),legend.background = element_rect(colour="grey80")) +
  scale_colour_grey(name = "Regime Assessment", 
                    labels = c("Government", "Human Rights","Democracy")) +
  scale_shape_discrete(name = "Model",
                       labels = c("Government", "Human Rights","Democracy"))


oneMod7 <- lm_robust(scale(govNonNA) ~ scale(SCI) + 
                      scale(gdppc) + scale(economicCondition) + scale(edu) + oil, 
                    data = fullData)
twoMod7 <- lm_robust(scale(hrNonNA) ~ scale(SCI) + 
                      scale(gdppc) + scale(economicCondition) + scale(edu) + oil, 
                    data = fullData)
threeMod7 <- lm_robust(scale(demNonNA) ~ scale(SCI) + 
                        scale(gdppc) + scale(economicCondition) + scale(edu) + oil, 
                      data = fullData)

#Figure 8

plot4NA <- dwplot(list(oneMod7,twoMod7,threeMod7),vline = geom_vline(xintercept = 0, colour = "grey60", linetype = 2)) %>%
  relabel_predictors(c("scale(SCI)"="SCI","scale(gdppc)"="GDPPC","scale(economicCondition)"="Economic Condition",
                       "scale(edu)"="Education","oil"="Oil")) +
  theme_bw() + ggtitle("Without FoE and Polity") + 
  theme(plot.title = element_text(hjust = 0.5),legend.background = element_rect(colour="grey80")) +
  scale_colour_grey(name = "Regime Assessment", 
                    labels = c("Government", "Human Rights","Democracy")) +
  scale_shape_discrete(name = "Model",
                       labels = c("Government", "Human Rights","Democracy"))

wvsReg1 <- ggarrange(plot1NA,plot2NA,plot3NA,plot4NA, common.legend = T, legend = "bottom")
wvsReg1 <- annotate_figure(wvsReg1,,
                           top = text_grob("OLS Regression Results (Nonresponse as NA)", 
                                           color = "black"))
ggexport(wvsReg1, filename = "Figure 8.pdf")

#Table A4

texreg(list(oneMod8,twoMod8,threeMod8,
          oneMod5,twoMod5,threeMod5,
          oneMod6,twoMod6,threeMod6,
          oneMod7,twoMod7,threeMod7),
       custom.coef.names = c("Intercept","SCI","Freedom of Expression","Polity","GDPPC", "Financial Situation", 
                             "Education","Oil"),
       custom.model.names  = c("Gov","HR","Dem","Gov","HR","Dem",
                               "Gov","HR","Dem","Gov","HR","Dem"),
       stars = c(0.05,0.01,0.001),
       reorder.coef = c(2,3,4,5,6,7,8,1),
       fontsize = "footnotesize",
       custom.header = list("FullSample"= 1:3,"FoE $>$ 0.5" = 4:6, "FoE $>$ Median" = 7:9, "w/o FoE and Polity" = 10:12),
       file = "TableA4.tex")



# Multilevel/Individual-level analysis




wvs <- merge(wvs,fullData,by = c("ccode","year"),all.x =T)

wvs$X003[wvs$X003==-1] <- NA
wvs$X003[wvs$X003==-2] <- NA
wvs$X003[wvs$X003==-3] <- NA
wvs$X003[wvs$X003==-4] <- NA
wvs$X003[wvs$X003==-5] <- NA
wvs$logAge <- log(wvs$X003)

wvs$X001[wvs$X001==-1] <- NA
wvs$X001[wvs$X001==-2] <- NA
wvs$X001[wvs$X001==-3] <- NA
wvs$X001[wvs$X001==-4] <- NA
wvs$X001[wvs$X001==-5] <- NA
wvs$Female <- wvs$X001 -1

wvs$X025[wvs$X025==-1] <- NA
wvs$X025[wvs$X025==-2] <- NA
wvs$X025[wvs$X025==-3] <- NA
wvs$X025[wvs$X025==-4] <- NA
wvs$X025[wvs$X025==-5] <- NA
wvs$Education <- (wvs$X025 > 1.5)*1 + (wvs$X025 > 5.5)*1

wvs$X049[wvs$X049==-1] <- NA
wvs$X049[wvs$X049==-2] <- NA
wvs$X049[wvs$X049==-3] <- NA
wvs$X049[wvs$X049==-4] <- NA
wvs$X049[wvs$X049==-5] <- NA
wvs$X049[wvs$X049==1] <- 0
wvs$X049[wvs$X049==2] <- 0
wvs$X049[wvs$X049==3] <- 1
wvs$X049[wvs$X049==4] <- 1
wvs$X049[wvs$X049==5] <- 1
wvs$X049[wvs$X049==6] <- 1
wvs$X049[wvs$X049==7] <- 1
wvs$X049[wvs$X049==8] <- 1
wvs$Urban <-  wvs$X049

wvs$A008[wvs$A008==1] <- 0
wvs$A008[wvs$A008==2] <- 1
wvs$A008[wvs$A008==-1] <- 2
wvs$A008[wvs$A008==-2] <- 2
wvs$A008[wvs$A008==-4] <- NA
wvs$A008[wvs$A008==-5] <- NA
wvs$A008 <- 4 - wvs$A008


wvs$logGDPPC <- log(as.numeric(wvs$gdppc))

wvs$polity2 <- wvs$polity2/20 +0.5

wvs$S021 <- as.factor(wvs$S021)

wvs$contributionSCI <- (wvs$govOB*wvs$hrOB*wvs$demOB*wvs$tvOB*wvs$lifeOB*wvs$trustOB)*(
  (wvs$govDK.x+wvs$demDK.x+wvs$hrDK.x)/3 - 
    (wvs$tvDK+wvs$lifeDK+wvs$trustDK)/3)


hmodel1 <- lmer(govMissNA ~ C006 + logAge + Female + Education + Urban + SCI + polity2 + logGDPPC + 
                  (1  | S021), data=wvs,na.action = na.omit)
hmodel2 <- lmer(hrMissNA ~ C006 + logAge + Female + Education + Urban + SCI + polity2 + logGDPPC + 
                  (1  | S021), data=wvs,na.action = na.omit)
hmodel3 <- lmer(demMissNA ~ C006 + logAge + Female + Education + Urban + SCI + polity2 + logGDPPC + 
                  (1  | S021), data=wvs,na.action = na.omit)
hmodel4 <- lmer(govMissNA ~ C006 + logAge + Female + Education + Urban + FoE + polity2 + logGDPPC + 
                  (1  | S021), data=wvs,na.action = na.omit)
hmodel5 <- lmer(hrMissNA ~ C006 + logAge + Female + Education + Urban + FoE + polity2 + logGDPPC + 
                  (1  | S021), data=wvs,na.action = na.omit)
hmodel6 <- lmer(demMissNA ~ C006 + logAge + Female + Education + Urban + FoE + polity2 + logGDPPC + 
                  (1  | S021), data=wvs,na.action = na.omit)



#Table A7

hmsum1 <- summary(hmodel1)
hmsum1
hmsum2 <- summary(hmodel2)
hmsum2
hmsum3 <- summary(hmodel3)
hmsum3
hmsum4 <- summary(hmodel4)
hmsum4
hmsum5 <- summary(hmodel5)
hmsum5
hmsum6 <- summary(hmodel6)
hmsum6


tableStar1 <- stargazer(hmodel1,hmodel2,hmodel3,hmodel4,hmodel5,hmodel6,
                        add.lines = list(c("Country-Waves",hmsum1$ngrps,hmsum2$ngrps,hmsum3$ngrps,
                                           hmsum4$ngrps,hmsum5$ngrps,hmsum6$ngrps)),
                        covariate.labels = c("Household Finances","log(Age)","Female","Education","Urban","SCI",
                                             "Freedom of Expansion","Polity","log(GDPPC)"),
                        dep.var.labels = c("Government Confidence","Human Rights","Democracy",
                                           "Government Confidence","Human Rights","Democracy"),
                        title = "LMER Results - SCI and Freedom of Expression (Manuscript)",
                        font.size = "footnotesize",
                        omit.stat = c("ll","aic","bic"), out = "TableA7.tex")

# Figure A11
plotNQ1 <- plot_ranef(hmodel1)
plotNQ2 <- plot_ranef(hmodel2)
plotNQ3 <- plot_ranef(hmodel3)
plotNQ4 <- plot_ranef(hmodel4)
plotNQ5 <- plot_ranef(hmodel5)
plotNQ6 <- plot_ranef(hmodel6)

diagnostics <- ggarrange(plotNQ1,plotNQ2,plotNQ3,plotNQ4,plotNQ5,plotNQ6, ncol = 2,nrow = 3,
                         labels = c("Government - SCI","Government - FoE","HR - SCI",
                         "HR - FoE","Democracy - SCI","Democracy - FoE"))

diagnostics

ggexport(diagnostics, filename = "FigureA11.pdf")

#Predicted Probabilities

predictSCI <- ggpredict(hmodel1, terms = c("SCI"))
predictSCIhr <- ggpredict(hmodel2, terms = c("SCI"))
predictSCIdem <- ggpredict(hmodel3, terms = c("SCI"))
predictFoE <- ggpredict(hmodel4, terms = c("FoE"))
predictFoEhr <- ggpredict(hmodel5, terms = c("FoE"))
predictFoEdem <- ggpredict(hmodel6, terms = c("FoE"))

this <- ggeffect(hmodel6, terms = c("FoE"))
plot(predictSCI, add.data=T)

# Plot the predictions 

SCIHM <- ggplot(predictSCI) + 
    geom_line(aes(x = x, y = predicted)) +          
    geom_ribbon(aes(x = x, ymin = predicted - 2*std.error, ymax = predicted + 2*std.error), 
                fill = "lightgrey", alpha = 0.5) + theme_bw() +
  ylab("Government") +
  xlab("") + ggtitle("SCI") + ylim(0,3)

SCIHMhr <- ggplot(predictSCIhr) + 
  geom_line(aes(x = x, y = predicted)) +          
  geom_ribbon(aes(x = x, ymin = predicted - 2*std.error, ymax = predicted + 2*std.error), 
              fill = "lightgrey", alpha = 0.5) + theme_bw() +
  ylab("Human Rights") +
  xlab("") + ylim(0,3)

SCIHMdem <- ggplot(predictSCIdem) + 
  geom_line(aes(x = x, y = predicted)) +         
  geom_ribbon(aes(x = x, ymin = predicted - 2*std.error, ymax = predicted + 2*std.error), 
              fill = "lightgrey", alpha = 0.5) + theme_bw() +
  ylab("Democracy") +
  xlab("") + ylim(0,10)


FoEHM <- ggplot(predictFoE) + 
  geom_line(aes(x = x, y = predicted)) +
  geom_ribbon(aes(x = x, ymin = predicted - 2*std.error, ymax = predicted + 2*std.error), 
              fill = "lightgrey", alpha = 0.6) + theme_bw() +
  ylab("") +
  xlab("")  + ggtitle("Freedom of Expression")+ ylim(0,3)

FoEHMhr <- ggplot(predictFoEhr) + 
  geom_line(aes(x = x, y = predicted)) + 
  geom_ribbon(aes(x = x, ymin = predicted - 2*std.error, ymax = predicted + 2*std.error), 
              fill = "lightgrey", alpha = 0.6) + theme_bw() +
  ylab("") +
  xlab("")  + ylim(0,3)

FoEHMhr

FoEHMdem <- ggplot(predictFoEdem) + 
  geom_line(aes(x = x, y = predicted)) +
  geom_ribbon(aes(x = x, ymin = predicted - 2*std.error, ymax = predicted + 2*std.error), 
              fill = "lightgrey", alpha = 0.6) + theme_bw() +
  ylab("") +
  xlab("")  + ylim(0,10)

FoEHMdem


#Table A9

hmodel7 <- lmer(govMissNA ~ C006 + logAge + Female + Education + Urban + SCI + 
                  (1  | S021), data=wvs,na.action = na.omit)
hmodel8 <- lmer(hrMissNA ~ C006 + logAge + Female + Education + Urban + SCI +  
                  (1  | S021), data=wvs,na.action = na.omit)
hmodel9 <- lmer(demMissNA ~ C006 + logAge + Female + Education + Urban + SCI +  
                  (1  | S021), data=wvs,na.action = na.omit)
hmodel10 <- lmer(govMissNA ~ C006 + logAge + Female + Education + Urban + FoE +  
                   (1  | S021), data=wvs,na.action = na.omit)
hmodel11<- lmer(hrMissNA ~ C006 + logAge + Female + Education + Urban + FoE + 
                  (1  | S021), data=wvs,na.action = na.omit)
hmodel12 <- lmer(demMissNA ~ C006 + logAge + Female + Education + Urban + FoE + 
                   (1  | S021), data=wvs,na.action = na.omit)


hmsum7 <- summary(hmodel7)
hmsum8 <- summary(hmodel8)
hmsum9 <- summary(hmodel9)
hmsum10 <- summary(hmodel10)
hmsum11 <- summary(hmodel11)
hmsum12 <- summary(hmodel12)

tableStar2 <- stargazer(hmodel7,hmodel8,hmodel9,hmodel10,hmodel11,hmodel12,
                        add.lines = list(c("Country-Waves",hmsum7$ngrps,hmsum8$ngrps,hmsum9$ngrps,
                                           hmsum10$ngrps,hmsum11$ngrps,hmsum12$ngrps)),
                        covariate.labels = c("Household Finances","log(Age)","Female","Education","Urban","SCI",
                                             "Freedom of Expansion"),
                        dep.var.labels = c("Government Confidence","Human Rights","Democracy",
                                           "Government Confidence","Human Rights","Democracy"),
                        title = "LMER Results - SCI and Freedom of Expression (Robustness Check)",
                        font.size = "footnotesize",
                        omit.stat = c("ll","aic","bic"), out = "TableA9.tex")


#Table A8

hmodel13 <- lmer(govMissNA ~ C006 + logAge + Female + Education + Urban + executive.comp + logGDPPC + 
                   (1  | S021), data=wvs,na.action = na.omit)
hmodel14 <- lmer(hrMissNA ~ C006 + logAge + Female + Education + Urban + executive.comp + logGDPPC + 
                   (1  | S021), data=wvs,na.action = na.omit)
hmodel15 <- lmer(demMissNA ~ C006 + logAge + Female + Education + Urban + executive.comp + logGDPPC +  
                   (1  | S021), data=wvs,na.action = na.omit)

hmsum13 <- summary(hmodel13)
hmsum14 <- summary(hmodel14)
hmsum15 <- summary(hmodel15)


stargazer(hmodel13,hmodel14,hmodel15,
          add.lines = list(c("Country-Waves",hmsum13$ngrps,hmsum14$ngrps,hmsum15$ngrps)),
          covariate.labels = c("Household Finances","log(Age)","Female","Education","Urban",
                               "Executive Competition", "log(GDPPC)"),
          dep.var.labels = c("Government Confidence","Human Rights","Democracy"),
          title = "LMER Results - Executive Competition (Manuscript)",
          omit.stat = c("ll","aic","bic"), out = "Table A8.tex")

predictEC <- ggpredict(hmodel13, terms = c("executive.comp [0:1]"))
predictEChr <- ggpredict(hmodel14, terms = c("executive.comp"))
predictECdem <- ggpredict(hmodel15, terms = c("executive.comp"))


#Executive Constraints Analysis


ECHM <- ggplot(predictEC,aes(x = x, y = predicted, ymin = conf.low, ymax = conf.high)) + 
  geom_line() +
  geom_errorbar(width=0.01) + geom_point() + 
  theme_bw()  + geom_line() +
  scale_x_continuous(breaks=c(0, 1)) +
  ylab("") +
  xlab("")  + ylim(0,3) +
  ggtitle("Executive Competition")
ECHM

ECHMhr <- ggplot(predictEChr,aes(x = x, y = predicted, ymin = conf.low, ymax = conf.high)) + 
  geom_line() +
  geom_errorbar(width=0.01) + geom_point() + 
  theme_bw()  + geom_line() +
  scale_x_continuous(breaks=c(0, 1)) +
  ylab(" ") +
  xlab(" ")  + ylim(0,3) +
  ggtitle(" ")
ECHMhr

ECHMdem <- ggplot(predictECdem,aes(x = x, y = predicted, ymin = conf.low, ymax = conf.high)) + 
  geom_line() +
  geom_errorbar(width=0.01) + geom_point() + 
  theme_bw()  + geom_line() +
  scale_x_continuous(breaks=c(0, 1)) +
  ylab("") +
  xlab("")  + ylim(1,10) +
  ggtitle("")
ECHMdem

#Figure A10

lmerplots <- ggarrange(SCIHM,FoEHM,ECHM,SCIHMhr,FoEHMhr,ECHMhr,SCIHMdem,FoEHMdem,ECHMdem, nrow=3,ncol=3)
ggexport(lmerplots,filename = "FigureA10.pdf")


#Table A10

hmodel30 <- lmer(E069_11 ~ C006 + logAge + Female + Education + Urban + SCI + polity2 + logGDPPC + 
                  (1  | S021), data=wvs,na.action = na.omit)
hmodel31 <- lmer(E124 ~ C006 + logAge + Female + Education + Urban + SCI + polity2 + logGDPPC + 
                  (1  | S021), data=wvs,na.action = na.omit)
hmodel32 <- lmer(E236 ~ C006 + logAge + Female + Education + Urban + SCI + polity2 + logGDPPC + 
                  (1  | S021), data=wvs,na.action = na.omit)
hmodel33 <- lmer(E069_11 ~ C006 + logAge + Female + Education + Urban + FoE + polity2 + logGDPPC + 
                  (1  | S021), data=wvs,na.action = na.omit)
hmodel34 <- lmer(E124 ~ C006 + logAge + Female + Education + Urban + FoE + polity2 + logGDPPC + 
                  (1  | S021), data=wvs,na.action = na.omit)
hmodel35 <- lmer(E236 ~ C006 + logAge + Female + Education + Urban + FoE + polity2 + logGDPPC + 
                  (1  | S021), data=wvs,na.action = na.omit)



hmsum30 <- summary(hmodel30)
hmsum31 <- summary(hmodel31)
hmsum32 <- summary(hmodel32)
hmsum33 <- summary(hmodel33)
hmsum34 <- summary(hmodel34)
hmsum35 <- summary(hmodel35)

tableStar3 <- stargazer(hmodel30,hmodel31,hmodel32,hmodel33,hmodel34,hmodel35,
                        add.lines = list(c("Country-Waves",hmsum1$ngrps,hmsum2$ngrps,hmsum3$ngrps,
                                           hmsum4$ngrps,hmsum5$ngrps,hmsum6$ngrps)),
                        covariate.labels = c("Household Finances","log(Age)","Female","Education","Urban","SCI",
                                             "Freedom of Expansion","Polity","log(GDPPC)"),
                        dep.var.labels = c("Government Confidence","Human Rights","Democracy",
                                           "Government Confidence","Human Rights","Democracy"),
                        title = "LMER Results - SCI and FoE (Robustness Check - Nonresponses coded as middle category)",
                        font.size = "footnotesize",
                        omit.stat = c("ll","aic","bic"), out = "TableA10.tex")


#Country Examinations
#Countries explored but not used in the analysis were left in for ease of exploration.
library(weights)



Uzbek <- data.frame("Country" = rep("Uzbekistan 2011",5), "Response"=c("None at all","Not very much", "Don't Know","Quite a lot", "A great deal"), 
                   "GovernmentSatisfaction" =wpct(wvs$E069_11[wvs$country=="Uzbekistan" & wvs$year==2011], 
                                                  wvs$S017,na.rm = T))
Uzbek
Iraq <- data.frame("Country" = rep("Iraq 2013",5), "Response"=c("None at all","Not very much", "Don't Know","Quite a lot", "A great deal"), 
                    "GovernmentSatisfaction" =wpct(wvs$E069_11[wvs$country=="Iraq" & wvs$year==2013], 
                                                   wvs$S017,na.rm = T))
Iraq
Mor11 <- data.frame("Country" = rep("Morocco 2011",5), "Response"=c("None at all","Not very much", "Don't Know","Quite a lot", "A great deal"), 
                    "GovernmentSatisfaction" =wpct(wvs$E069_11[wvs$country=="Morocco" & wvs$year==2011], 
                                                   wvs$S017,na.rm = T))
Mor11

Iran <- data.frame("Country" = rep("Iran 2007",5), "Response"=c("None at all","Not very much", "Don't Know","Quite a lot", "A great deal"), 
                    "GovernmentSatisfaction" =wpct(wvs$E069_11[wvs$country=="Iran" & wvs$year==2007], 
                                                   wvs$S017,na.rm = T))
Iran
Azer <- data.frame("Country" = rep("Azerbaijan 2011",5), "Response"=c("Don't Know","None at all","Not very much", "Quite a lot", "A great deal"), 
                   "GovernmentSatisfaction" =c(0,wpct(wvs$E069_11[wvs$country=="Azerbaijan" & wvs$year==2011], 
                                                  wvs$S017,na.rm = T)))
Azer
France <- data.frame("Country" = rep("France 2006",5), "Response"=c("None at all","Not very much", "Don't Know","Quite a lot", "A great deal"), 
                   "GovernmentSatisfaction" =wpct(wvs$E069_11[wvs$country=="France" & wvs$year==2006], 
                                                  wvs$S017,na.rm = T))
France
UK <- data.frame("Country" = rep("UK 2005",5), "Response"=c("None at all","Not very much", "Don't Know","Quite a lot", "A great deal"), 
                     "GovernmentSatisfaction" =wpct(wvs$E069_11[wvs$country=="United Kingdom" & wvs$year==2005], 
                                                    wvs$S017,na.rm = T))
UK
Pak <- data.frame("Country" = rep("Pakistan 2001",5), "Response"=c("None at all","Not very much", "Don't Know","Quite a lot", "A great deal"), 
                 "GovernmentSatisfaction" =wpct(wvs$E069_11[wvs$country=="Pakistan" & wvs$year==2001], 
                                                wvs$S017,na.rm = T))
Pak
Kyr <- data.frame("Country" = rep("Kyrgyzstan 2003",5), "Response"=c("None at all","Not very much", "Don't Know","Quite a lot", "A great deal"), 
                  "GovernmentSatisfaction" =wpct(wvs$E069_11[wvs$country=="Kyrgyzstan" & wvs$year==2003], 
                                                 wvs$S017,na.rm = T))
Kyr



CountryMerge <- rbind(UK, Iraq,Uzbek, Pak)
CountryMerge

CountryMerge$Response <- ordered(CountryMerge$Response, levels=c("A great deal","Quite a lot","Not very much","None at all","Don't Know"))
CountryMerge$Country <- ordered(CountryMerge$Country, levels=c("UK 2005","Iraq 2013","Uzbekistan 2011","Pakistan 2001"))


fullData$FoE[fullData$country=="Pakistan" & fullData$year==2001]
fullData$FoE[fullData$country=="United Kingdom" & fullData$year==2005]
fullData$FoE[fullData$country=="Uzbekistan" & fullData$year==2011]
fullData$FoE[fullData$country=="Iraq" & fullData$year==2013]

#Figure 6

countryPlot <- ggplot(CountryMerge,aes(x=Country, y=100*GovernmentSatisfaction,fill=Response)) + 
  geom_bar(stat = "identity") + ylab("Percentage") + xlab("") + 
  ggtitle("Confidence in Government") + theme(legend.title = element_blank()) +
  geom_segment(aes(x=0.55,y=105.8,xend=3.45,yend=105.8),col="dark green", size=1.2) +
  geom_segment(aes(x=3.55,y=105.8,xend=4.45,yend=105.8),col="red", size=1.2) +
  annotate("text",x=2,y=108.7,col="dark green", label = "Low SCI") +
  annotate("text",x=4,y=108.7,col="red", label = "High SCI") +
  annotate("text",x=0.7,y=102.2,col="black", label = "FoE:") +
  annotate("text",x=1,y=102.2,col="black", label = "0.96") +
  annotate("text",x=2,y=102.2,col="black", label = "0.67") +
  annotate("text",x=3,y=102.2,col="black", label = "0.04") +
  annotate("text",x=4,y=102.2,col="black", label = "0.77") +
  scale_y_continuous(breaks=seq(0,100,by=25))
  
  
countryPlot

ggexport(countryPlot, filename="Figure6.pdf")


#TableA3
summarySCI <- data.frame("Country" = fullData$Country,"Year" = as.character(fullData$year),"SCI" = fullData$SCI,"SCIall" = fullData$SCIallPos)
summarySCI <- summarySCI[order(summarySCI$Country, summarySCI$Year),]

summarySCI <- summarySCI[!is.na(summarySCI$SCIall),]

stargazer(summarySCI,summary = F,rownames = F, 
          covariate.labels = c("Country","Year","SCI - Complete Data","SCI - Any Available"),
          title="SCI Scores for the Empirical Analysis", out = "TableA3.tex")

#Table A2
summaryData <- data.frame("Gov"=fullData$govNonNA,"HR"=fullData$hrNonNA,
                          "Dem"=fullData$demNonNA,"SCI"=fullData$SCI,
                          "FoE"=fullData$FoE,
                          "Polity"=fullData$polity2,
                          "GDPPC"=fullData$gdppc/1000,
                          "Economic"=fullData$economicCondition,
                          "Education"=fullData$edu,"Oil"=fullData$oil)

stargazer(summaryData,
          covariate.labels = c("Government Confidence","Human Rights","Democracy (Survey)","SCI",
                               "Freedom of Expression","Polity","GDPPC (thousands)","Household Financial Situation","Education",
                               "Oil"),
          title="Summary Statistics (Aggregate Data)",
          font.size = "small",
          out = "TableA2.tex")


#OLS with clustered standard errors and no country-year fixed effects

IndClus1 <- lm_robust(govMissNA ~ C006 + logAge + Female + Education + Urban + SCI + polity2 + logGDPPC,
                     clusters = S021, se_type = "CR2",  data=wvs)
IndClus2 <- lm_robust(hrMissNA ~ C006 + logAge + Female + Education + Urban + SCI + polity2 + logGDPPC,
                      clusters = S021, se_type = "CR2",  data=wvs)
IndClus3 <- lm_robust(demMissNA ~ C006 + logAge + Female + Education + Urban + SCI + polity2 + logGDPPC,
                      clusters = S021, se_type = "CR2",  data=wvs)
IndClus4 <- lm_robust(govMissNA ~ C006 + logAge + Female + Education + Urban + FoE + polity2 + logGDPPC,
                      clusters = S021, se_type = "CR2",  data=wvs)
IndClus5 <- lm_robust(hrMissNA ~ C006 + logAge + Female + Education + Urban + FoE + polity2 + logGDPPC,
                      clusters = S021, se_type = "CR2",  data=wvs)
IndClus6 <- lm_robust(demMissNA ~ C006 + logAge + Female + Education + Urban + FoE + polity2 + logGDPPC,
                      clusters = S021, se_type = "CR2",  data=wvs)
IndClus7 <- lm_robust(govMissNA ~ C006 + logAge + Female + Education + Urban + executive.comp + polity2 + logGDPPC,
                      clusters = S021, se_type = "CR2",  data=wvs)
IndClus8 <- lm_robust(hrMissNA ~ C006 + logAge + Female + Education + Urban + executive.comp + polity2 + logGDPPC,
                      clusters = S021, se_type = "CR2",  data=wvs)
IndClus9 <- lm_robust(demMissNA ~ C006 + logAge + Female + Education + Urban + executive.comp + polity2 + logGDPPC,
                      clusters = S021, se_type = "CR2",  data=wvs)


predictSCIic <- ggpredict(IndClus1, terms = c("SCI"))
predictSCIhric <- ggpredict(IndClus2, terms = c("SCI"))
predictSCIdemic <- ggpredict(IndClus3, terms = c("SCI"))
predictFoEic <- ggpredict(IndClus4, terms = c("FoE"))
predictFoEhric <- ggpredict(IndClus5, terms = c("FoE"))
predictFoEdemic <- ggpredict(IndClus6, terms = c("FoE"))
predictECic <- ggpredict(IndClus7, terms = c("executive.comp"))
predictEChric <- ggpredict(IndClus8, terms = c("executive.comp"))
predictECdemic <- ggpredict(IndClus9, terms = c("executive.comp"))


SCIOLSgov <- ggplot(predictSCIic) + 
  geom_line(aes(x = x, y = predicted)) +        
  geom_ribbon(aes(x = x, ymin = conf.low, ymax = conf.high), 
              fill = "lightgrey", alpha = 0.5) + theme_bw() +
  ylab("Government") +
  xlab("") + ggtitle("SCI") + ylim(0,3)

SCIOLShr <- ggplot(predictSCIhric) + 
  geom_line(aes(x = x, y = predicted)) + 
  geom_ribbon(aes(x = x, ymin = conf.low, ymax = conf.high), 
              fill = "lightgrey", alpha = 0.5) + theme_bw() +
  ylab("Human Rights") +
  xlab("") + ylim(0,3)

SCIOLSdem <- ggplot(predictSCIdemic) + 
  geom_line(aes(x = x, y = predicted)) +  
  geom_ribbon(aes(x = x, ymin = conf.low, ymax = conf.high), 
              fill = "lightgrey", alpha = 0.5) + theme_bw() +
  ylab("Democracy") +
  xlab("") + ylim(0,10)


FoEOLSgov <- ggplot(predictFoEic) + 
  geom_line(aes(x = x, y = predicted)) +
  geom_ribbon(aes(x = x, ymin = conf.low, ymax = conf.high), 
              fill = "lightgrey", alpha = 0.6) + theme_bw() +
  ylab("") +
  xlab("")  + ggtitle("Freedom of Expression")+ ylim(0,3)

FoEOLShr <- ggplot(predictFoEhric) + 
  geom_line(aes(x = x, y = predicted)) + 
  geom_ribbon(aes(x = x, ymin = conf.low, ymax = conf.high), 
              fill = "lightgrey", alpha = 0.6) + theme_bw() +
  ylab("") +
  xlab("")  + ylim(0,3)

FoEOLShr

FoEOLSdem <- ggplot(predictFoEdemic) + 
  geom_line(aes(x = x, y = predicted)) +
  geom_ribbon(aes(x = x, ymin = conf.low, ymax = conf.high), 
              fill = "lightgrey", alpha = 0.6) + theme_bw() +
  ylab("") +
  xlab("")  + ylim(0,10)

FoEOLSdem

ECOLSgov <- ggplot(predictECic,aes(x = x, y = predicted, ymin = conf.low, ymax = conf.high)) + 
  geom_line() +
  geom_errorbar(width=0.01) + geom_point() + 
  theme_bw()  + geom_line() +
  scale_x_continuous(breaks=c(0, 1)) +
  ylab("") +
  xlab("")  + ylim(0,3) +
  ggtitle("Executive Competition")


ECOLShr <- ggplot(predictEChric,aes(x = x, y = predicted, ymin = conf.low, ymax = conf.high)) + 
  geom_line() +
  geom_errorbar(width=0.01) + geom_point() + 
  theme_bw()  + geom_line() +
  scale_x_continuous(breaks=c(0, 1)) +
  ylab("") +
  xlab("")  + ylim(0,3) +
  ggtitle("")

ECOLSdem <- ggplot(predictECdemic,aes(x = x, y = predicted, ymin = conf.low, ymax = conf.high)) + 
  geom_line() +
  geom_errorbar(width=0.01) + geom_point() + 
  theme_bw()  + geom_line() +
  scale_x_continuous(breaks=c(0, 1)) +
  ylab("") +
  xlab("")  + ylim(0,10) +
  ggtitle("")


# Figure A9

figsOLS <- ggarrange(SCIOLSgov,FoEOLSgov,ECOLSgov,
                     SCIOLShr,FoEOLShr,ECOLShr,
                     SCIOLSdem,FoEOLSdem,ECOLSdem,
                     ncol=3,nrow=3)
ggexport(figsOLS,dpi=900,filename = "FigureA9.pdf")

#PC Regression

wvs1 <- wvs %>% select(govMissNA, C006, logAge, Female, Education, Urban, SCI, polity2, logGDPPC, S021)
wvs1 <- wvs1[complete.cases(wvs1),]
wvs2 <- wvs %>% select(hrMissNA, C006, logAge, Female, Education, Urban, SCI, polity2, logGDPPC, S021)
wvs2 <- wvs2[complete.cases(wvs2),]
wvs3 <- wvs %>% select(demMissNA, C006, logAge, Female, Education, Urban, SCI, polity2, logGDPPC, S021)
wvs3 <- wvs3[complete.cases(wvs3),]
wvs4 <- wvs %>% select(govMissNA, C006, logAge, Female, Education, Urban, FoE, polity2, logGDPPC, S021)
wvs4 <- wvs4[complete.cases(wvs4),]
wvs5 <- wvs %>% select(hrMissNA, C006, logAge, Female, Education, Urban, FoE, polity2, logGDPPC, S021)
wvs5 <- wvs5[complete.cases(wvs5),]
wvs6 <- wvs %>% select(demMissNA, C006, logAge, Female, Education, Urban, FoE, polity2, logGDPPC, S021)
wvs6 <- wvs6[complete.cases(wvs6),]
wvs7 <- wvs %>% select(govMissNA, C006, logAge, Female, Education, Urban, executive.comp, logGDPPC, S021)
wvs7 <- wvs7[complete.cases(wvs7),]
wvs8 <- wvs %>% select(hrMissNA, C006, logAge, Female, Education, Urban, executive.comp, logGDPPC, S021)
wvs8 <- wvs8[complete.cases(wvs8),]
wvs9 <- wvs %>% select(demMissNA, C006, logAge, Female, Education, Urban, executive.comp, logGDPPC, S021)
wvs9 <- wvs9[complete.cases(wvs9),]



FEReg1 <- lm_robust(govMissNA ~ C006 + logAge + Female + Education + Urban, fixed_effects = S021,
                    clusters = S021, se_type = "CR2",  data=wvs1)
FEReg2 <- lm_robust(hrMissNA ~ C006 + logAge + Female + Education + Urban, fixed_effects = S021,
                    clusters = S021, se_type = "CR2",  data=wvs2)
FEReg3 <- lm_robust(demMissNA ~ C006 + logAge + Female + Education + Urban, fixed_effects = S021,
                    clusters = S021, se_type = "CR2",  data=wvs3)
FEReg4 <- lm_robust(govMissNA ~ C006 + logAge + Female + Education + Urban, fixed_effects = S021,
                    clusters = S021, se_type = "CR2",  data=wvs4)
FEReg5 <- lm_robust(hrMissNA ~ C006 + logAge + Female + Education + Urban, fixed_effects = S021,
                    clusters = S021, se_type = "CR2",  data=wvs5)
FEReg6 <- lm_robust(demMissNA ~ C006 + logAge + Female + Education + Urban, fixed_effects = S021,
                    clusters = S021, se_type = "CR2",  data=wvs6)
FEReg7 <- lm_robust(govMissNA ~ C006 + logAge + Female + Education + Urban, fixed_effects = S021,
                    clusters = S021, se_type = "CR2",  data=wvs7)
FEReg8 <- lm_robust(hrMissNA ~ C006 + logAge + Female + Education + Urban, fixed_effects = S021,
                    clusters = S021, se_type = "CR2",  data=wvs8)
FEReg9 <- lm_robust(demMissNA ~ C006 + logAge + Female + Education + Urban, fixed_effects = S021,
                    clusters = S021, se_type = "CR2",  data=wvs9)


forCoef1 <- lm(govMissNA ~ C006 + logAge + Female + Education + Urban + factor(S021), data=wvs1)
forCoef2 <- lm(hrMissNA ~ C006 + logAge + Female + Education + Urban + factor(S021), data=wvs2)
forCoef3 <- lm(demMissNA ~ C006 + logAge + Female + Education + Urban + factor(S021), data=wvs3)
forCoef4 <- lm(govMissNA ~ C006 + logAge + Female + Education + Urban + factor(S021), data=wvs4)
forCoef5 <- lm(hrMissNA ~ C006 + logAge + Female + Education + Urban + factor(S021), data=wvs5)
forCoef6 <- lm(demMissNA ~ C006 + logAge + Female + Education + Urban + factor(S021), data=wvs6)
forCoef7 <- lm(govMissNA ~ C006 + logAge + Female + Education + Urban + factor(S021), data=wvs7)
forCoef8 <- lm(hrMissNA ~ C006 + logAge + Female + Education + Urban + factor(S021), data=wvs8)
forCoef9 <- lm(demMissNA ~ C006 + logAge + Female + Education + Urban + factor(S021), data=wvs9)

#Table A6

stargazer(forCoef1, forCoef2, forCoef3, forCoef4, forCoef5, forCoef6, forCoef7, forCoef8, forCoef9, 
          se = starprep(FEReg1, FEReg2, FEReg3, FEReg4, FEReg5, FEReg6, FEReg7, FEReg8, FEReg9),
          keep = c("C006", "log(Age)", "Female", "Education", "Urban"),
          covariate.labels = c("Household Finances","log(Age)","Female","Education","Urban"),
          dep.var.labels = c("GC","HR","Dem",
                             "GC","HR","Dem",
                             "GC","HR","Dem"),
          title = "First Step of Per Cluster Regression Analysis",
          star.cutoffs = NA,
          font.size = "scriptsize",
          omit.stat = c("f","adj.rsq","rsq","ser"),
          notes.append = F,
          notes = "(Standard errors in parentheses)",
          out = "Table A6.tex")

library(texreg)
texreg(list(FEReg1, FEReg2, FEReg3), include.ci=F)

wvs$govFit1 <- wvs$govMissNA - wvs$C006*FEReg1$coefficients[1] - wvs$logAge*FEReg1$coefficients[2] - 
  wvs$Female*FEReg1$coefficients[3] - wvs$Education*FEReg1$coefficients[4] - wvs$Urban*FEReg1$coefficients[5]
wvs$hrFit1 <- wvs$hrMissNA - wvs$C006*FEReg2$coefficients[1] - wvs$logAge*FEReg2$coefficients[2] - 
  wvs$Female*FEReg2$coefficients[3] - wvs$Education*FEReg2$coefficients[4] - wvs$Urban*FEReg2$coefficients[5]
wvs$demFit1 <- wvs$demMissNA - wvs$C006*FEReg3$coefficients[1] - wvs$logAge*FEReg3$coefficients[2] - 
  wvs$Female*FEReg3$coefficients[3] - wvs$Education*FEReg3$coefficients[4] - wvs$Urban*FEReg3$coefficients[5]
wvs$govFit2 <- wvs$govMissNA - wvs$C006*FEReg4$coefficients[1] - wvs$logAge*FEReg4$coefficients[2] - 
  wvs$Female*FEReg4$coefficients[3] - wvs$Education*FEReg4$coefficients[4] - wvs$Urban*FEReg4$coefficients[5]
wvs$hrFit2 <- wvs$hrMissNA - wvs$C006*FEReg5$coefficients[1] - wvs$logAge*FEReg5$coefficients[2] - 
  wvs$Female*FEReg5$coefficients[3] - wvs$Education*FEReg5$coefficients[4] - wvs$Urban*FEReg5$coefficients[5]
wvs$demFit2 <- wvs$demMissNA - wvs$C006*FEReg6$coefficients[1] - wvs$logAge*FEReg6$coefficients[2] - 
  wvs$Female*FEReg6$coefficients[3] - wvs$Education*FEReg6$coefficients[4] - wvs$Urban*FEReg6$coefficients[5]
wvs$govFit3 <- wvs$govMissNA - wvs$C006*FEReg7$coefficients[1] - wvs$logAge*FEReg7$coefficients[2] - 
  wvs$Female*FEReg7$coefficients[3] - wvs$Education*FEReg7$coefficients[4] - wvs$Urban*FEReg7$coefficients[5]
wvs$hrFit3 <- wvs$hrMissNA - wvs$C006*FEReg8$coefficients[1] - wvs$logAge*FEReg8$coefficients[2] - 
  wvs$Female*FEReg8$coefficients[3] - wvs$Education*FEReg8$coefficients[4] - wvs$Urban*FEReg8$coefficients[5]
wvs$demFit3 <- wvs$demMissNA - wvs$C006*FEReg9$coefficients[1] - wvs$logAge*FEReg9$coefficients[2] - 
  wvs$Female*FEReg9$coefficients[3] - wvs$Education*FEReg9$coefficients[4] - wvs$Urban*FEReg9$coefficients[5]

feReg <- wvs %>%
  group_by(S021,ccode,year) %>%
  summarise_at(c("govFit1","demFit1","hrFit1", "govFit2","demFit2","hrFit2", "govFit3","demFit3","hrFit3", "SCI", "polity2", "logGDPPC", "FoE", "executive.comp"), mean,na.rm=T)

finalGov <- lm_robust(govFit1 ~ SCI + polity2 + logGDPPC, data = feReg)
summary(finalGov)
finalDem <- lm_robust(demFit1 ~ SCI + polity2 + logGDPPC, data = feReg)
summary(finalDem)
finalHr <- lm_robust(hrFit1 ~ SCI + polity2 + logGDPPC, data = feReg)
summary(finalHr)
finalGovfe <- lm_robust(govFit2 ~ FoE + polity2 + logGDPPC, data = feReg)
summary(finalGovfe)
finalDemfe <- lm_robust(demFit2 ~ FoE + polity2 + logGDPPC, data = feReg)
summary(finalDemfe)
finalHrfe <- lm_robust(hrFit2 ~ FoE + polity2 + logGDPPC, data = feReg)
summary(finalHrfe)
finalGovec <- lm_robust(govFit3 ~ executive.comp + logGDPPC, data = feReg)
summary(finalGovec)
finalDemec <- lm_robust(demFit3 ~ executive.comp + logGDPPC, data = feReg)
summary(finalDemec)
finalHrec <- lm_robust(hrFit3 ~ executive.comp + logGDPPC, data = feReg)
summary(finalHrec)

predictSCIfe <- ggeffect(finalGov, terms = c("SCI"))
predictSCIhrfe <- ggeffect(finalHr, terms = c("SCI"))
predictSCIdemfe <- ggeffect(finalDem, terms = c("SCI"))
predictFoEfe <- ggeffect(finalGovfe, terms = c("FoE"))
predictFoEhrfe <- ggeffect(finalHrfe, terms = c("FoE"))
predictFoEdemfe <- ggeffect(finalDemfe, terms = c("FoE"))
predictECfe <- ggeffect(finalGovec, terms = c("executive.comp"))
predictEChrfe <- ggeffect(finalHrec, terms = c("executive.comp"))
predictECdemfe <- ggeffect(finalDemec, terms = c("executive.comp"))


SCIFEgov <- ggplot(predictSCIfe) + 
  geom_line(aes(x = x, y = predicted)) +
  geom_ribbon(aes(x = x, ymin = conf.low, ymax = conf.high), 
              fill = "lightgrey", alpha = 0.5) + theme_bw() +
  ylab("Government") +
  xlab("") + ggtitle("SCI") + ylim(0,3)

SCIFEhr <- ggplot(predictSCIhrfe) + 
  geom_line(aes(x = x, y = predicted)) + 
  geom_ribbon(aes(x = x, ymin = conf.low, ymax = conf.high), 
              fill = "lightgrey", alpha = 0.5) + theme_bw() +
  ylab("Human Rights") +
  xlab("") + ylim(0,3)

SCIFEdem <- ggplot(predictSCIdemfe) + 
  geom_line(aes(x = x, y = predicted)) +
  geom_ribbon(aes(x = x, ymin = conf.low, ymax = conf.high), 
              fill = "lightgrey", alpha = 0.5) + theme_bw() +
  ylab("Democracy") +
  xlab("") + ylim(0,10)


FoEFEgov <- ggplot(predictFoEfe) + 
  geom_line(aes(x = x, y = predicted)) + 
  geom_ribbon(aes(x = x, ymin = conf.low, ymax = conf.high), 
              fill = "lightgrey", alpha = 0.6) + theme_bw() +
  ylab("") +
  xlab("")  + ggtitle("Freedom of Expression")+ ylim(0,3)

FoEFEhr <- ggplot(predictFoEhrfe) + 
  geom_line(aes(x = x, y = predicted)) + 
  geom_ribbon(aes(x = x, ymin = conf.low, ymax = conf.high), 
              fill = "lightgrey", alpha = 0.6) + theme_bw() +
  ylab("") +
  xlab("")  + ylim(0,3)

FoEFEhr

FoEFEdem <- ggplot(predictFoEdemfe) + 
  geom_line(aes(x = x, y = predicted)) + 
  geom_ribbon(aes(x = x, ymin = conf.low, ymax = conf.high), 
              fill = "lightgrey", alpha = 0.6) + theme_bw() +
  ylab("") +
  xlab("")  + ylim(0,10)

FoEFEdem

ECFEgov <- ggplot(predictECfe,aes(x = x, y = predicted, ymin = conf.low, ymax = conf.high)) + 
  geom_line() +
  geom_errorbar(width=0.01) + geom_point() + 
  theme_bw()  + geom_line() +
  scale_x_continuous(breaks=c(0, 1)) +
  ylab("") +
  xlab("")  + ylim(0,3) +
  ggtitle("Executive Competition")


ECFEhr <- ggplot(predictEChrfe,aes(x = x, y = predicted, ymin = conf.low, ymax = conf.high)) + 
  geom_line() +
  geom_errorbar(width=0.01) + geom_point() + 
  theme_bw()  + geom_line() +
  scale_x_continuous(breaks=c(0, 1)) +
  ylab("") +
  xlab("")  + ylim(0,3) +
  ggtitle("")

ECFEdem <- ggplot(predictECdemfe,aes(x = x, y = predicted, ymin = conf.low, ymax = conf.high)) + 
  geom_line() +
  geom_errorbar(width=0.01) + geom_point() + 
  theme_bw()  + geom_line() +
  scale_x_continuous(breaks=c(0, 1)) +
  ylab("") +
  xlab("")  + ylim(0,10) +
  ggtitle("")

#Figure 9

figsFE <- ggarrange(SCIFEgov,FoEFEgov,ECFEgov,
                     SCIFEhr,FoEFEhr,ECFEhr,
                     SCIFEdem,FoEFEdem,ECFEdem,
                     ncol=3,nrow=3)
ggexport(figsFE,dpi=900,filename = "Figure9.pdf")

